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INVISCID RADIATING SHOCK LAYERS ABOUT SPHERES 


TRAVELING AT HYPERBOLIC SPEEDS IN AIR* 

By Richard W. Barnwell 
Langley Research Center 

SUMMARY 

Time-dependent finite-difference techniques are used to obtain numerical solutions 
for the problem of the inviscid flow of radiating equilibrium air past spheres traveling at 
hyperbolic speeds. The effects of absorption are included, and results are presented for 
both gray and nongray absorption coefficient models for spheres with different radii. It 
is shown that the nondimensional heat-flux distributions for the gray and nongray models 
are similar and that these distributions are weak functions of the radius of the sphere and 
the altitude and strong functions of the flight velocity. 

INTRODUCTION 

The problem of radiating flow past blunt bodies traveling at hyperbolic speeds in air 
has been the subject of intensive theoretical investigation during the past decade. Solu- 
tions have been obtained for inviscid and viscous flows and for both transparent and 
absorbing radiation models. In reference 1, Goulard has investigated transparent, invis- 
cid flow fields. Thomas (ref. 2) has extended Goulard's work to include the effects of 
absorption for cases for which the absorption coefficient is small by replacing the absorp- 
tion integral at a point in the flow field with a Taylor series expansion about the point. An 
iterative numerical procedure is used by Howe and Viegas (ref. 3) to treat viscous flow 
fields which emit and absorb thermal radiation. Olstad (refs. 4 and 5) employs the 
Poincare- Lighthill-Kuo method to obtain perturbation solutions to the inviscid radiation 
problem. Callis (ref. 6) solves this same problem with a time-dependent finite-difference 
technique of second-order accuracy. In all the references cited, modified one-dimensional 
flow fields are considered, and the solutions which are obtained are applicable in the 
vicinity of the stagnation streamline. 


*An earlier version of part of the analysis which is presented in this paper was used 
as a thesis in partial fulfillment of the requirements for the degree of Doctor of Philosophy 
in Engineering Mechanics, Virginia Polytechnic Institute, Blacksburg, Virginia, June 1968. 


Solutions which are not restricted to the stagnation streamline are presented in 
references 7 to 10. In references 7 and 8, Hoshizaki and Wilson use an integral method 
to obtain solutions for inviscid and viscous flows, respectively, for a transparent radia- 
tion model. These authors extend their work to include the effects of absorption in ref- 
erence 9. Cheng and Vincenti (ref. 10) use an inverse method to treat the problem of 
inviscid flow past blunt bodies and consider the effects of absorption. 

The frequency dependence of the absorption coefficient for high-temperature air is 
considered in references 4, 5, 6, 9, and 11. Hoshizaki and Wilson (ref. 9) perform the 
integration with respect to frequency in the radiative heat-transfer term of the energy 
equation numerically, whereas Olstad (refs. 4 and 5) and Smith and Hassan (ref. 11) use 
step functions to approximate the frequency-dependent absorption coefficient and perform 
the integration with respect to frequency analytically. The step-function absorption coef- 
ficient models are constructed by dividing the spectrum into intervals and using mean 
values of the absorption coefficients for each interval. 

In the present paper, the computation of the radiating shock layer about a sphere is 
made by using a time-dependent method of first-order accuracy. This method is a "dis- 
crete shock" method which is similar to that of Moretti and Abbett (ref. 12). The shock 
wave is treated as a boundary at which the flow properties are discontinuous. With time- 
dependent, finite-difference methods, one can avoid the assumption of functional forms 
for some of the normal flow property profiles as is done with the integral method which 
is used in references 7, 8, and 9. Another advantage is that time-dependent methods can 
be used to solve the direct blunt-body radiation problem as opposed to the inverse prob- 
lem, which is solved in reference 10. 

The problem of time-dependent supersonic flow about a blunt body is treated as an 
initial-value problem. By starting with an initial approximate solution, the solution is 
determined at successive time steps until the steady solution is approached asymptoti- 
cally. The inclusion of the partial derivatives with respect to time in the equations for 
fluid flow makes those equations hyperbolic at all points in the flow field including regions 
where the flow is subsonic. This feature makes the time-dependent approach an attrac- 
tive one for computing supersonic flow past blunt bodies. The governing partial differen- 
tial equations are approximated with finite-difference equations which are solved with 
forward- marching techniques. 

The time -dependent approach is particularly applicable to the problem of radiating 
flow with absorption because it provides a straightforward means of solving the energy 
equation, which is an integro-differential equation. When determining the solution at a 
new time step, the known solution at the previous step is used to evaluate the integral 
term. 
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As was indicated earlier, the present paper treats the inviscid flow fields about 
spheres traveling at hyperbolic velocities at high altitudes in the earth's atmosphere. 
Spheres with several different radii are treated, and both frequency-dependent and 
frequency-independent absorption coefficient models are used. The frequency-dependent 
absorption coefficient is approximated with a step-function model which is similar to 
those used in references 4, 5, 6, and 11. It is assumed that the gas in the free stream 
neither emits nor absorbs thermal energy and that the body surface neither emits nor 
reflects. For the purpose of calculating the absorption integral in the energy equation, 
the radiating shock layer is treated as an infinite slab with the same width as the local- 
shock-layer thickness. The thermodynamic state of the gas is described by curve fits of 
previously calculated equilibrium air data. 

SYMBOLS 

Ai,Bi,Ci,Di,Ei,Fi quantities defined by equations (10) 

An[j,s] integral defined by equation (E9) 

Aj/[r,s] optical depth defined by equation (40) 
a speed of sound 

a Q reference speed of sound, 331.8 m/sec 

aN,b N exponents expressing power -law variation of absorption coefficient with 
density and temperature, respectively 

B^f frequency integral of Planck function defined by equation (E8) 

Bj, Planck function 

C constant used in equation (37) 

C quantity defined by equation (E26) 

c speed of light 

d exponent expressing power- law variation of shock- layer thickness with cos 0 
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E 


total energy per unit mass 


E n (z) exponential integral defined by equation (39) 

e internal energy per unit mass 


F(z) function defined by equation (43) 



hypergeometric series defined by equation (E28) 


f scale constant used in equation (E16) 

f constant defined by equation (E18) 

g exponent expressing power-law variation of Planck function with hc/kXT 

H total enthalpy per unit mass 

h static enthalpy per unit mass, also Planck constant 

h^ scale constant for 77 coordinate 

h altitude 


Ij, radiation intensity 

1 ^, 1 ^, radiation intensity away from body surface and toward it, respectively 

K shock curvature 


k Boltzmann constant 


Z spatial position, r = r^ + (Z - l)Ar and Y = (Z - 1)AY 

m spatial position, 0 = (m - 1)A0 and 0 = (m - 1)A0 

M,N,n integers 

/V 

N quantity defined by equation (E25) 
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p pressure 

p Q reference pressure, 1.013 x 10^ newtons/m^ 

Q total radiative heat input per unit volume and time defined by equation (D15) 

Q^j- quantity defined by equation (E5) 

monochromatic radiative heat input per unit volume and time defined by 
equation (D9) 

q total radiative heat flux defined by equation (D8) 

q^ monochromatic radiative heat flux defined by equation (D3) 

R perpendicular distance from axis of symmetry to a point in flow field, also 

gas constant for air, 2.870 x 10^ joules/kg-°K 

R 0 quantity defined by equation (B7) 

R,S exponents expressing power-law variation of temperature with pressure and 

density, respectively 

r radial coordinate 

S entropy 

T temperature 

T Q reference temperature, 273.2° K 

t time 

U,V velocity components normal and parallel to shock wave, respectively 

u,v velocity components normal and parallel to body surface, respectively 

magnitude of free-stream velocity 
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component of shock velocity normal to body surface defined by equation (14) 
Cartesian coordinate, also Goulard coordinate 

nondimensional coordinate normal to body surface defined by equation (11c) 
coordinate normal to body surface, also Cartesian coordinate 
nondimensional coordinate defined by equation (E7) 
length along line of propagation of radiation 
absorption coefficient 

mean absorption coefficient for interval ^ v ^ ^N+l 
Planck mean absorption coefficient 

acute angle between normal to shock wave and free-stream direction 

energy depletion parameter defined by equation (46) 

energy depletion parameter defined by equation (E22) 

ratio of specific heats 

ratio of static enthalpy to internal energy 

mesh spacings for time t and coordinates r and 0, respectively 
mesh spacings for time r and coordinates Y and ©, respectively 

shock-layer thickness along a line normal to body surface defined by 
equation (12) 

shock-detachment distance for adiabatic flow 

acute angle between normals to shock wave and body surface 
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£ quantity defined by equation (E23) 

77 coordinate along shock wave 

0,0 polar angle 

X wavelength 

l-L quantity defined by equations (9) 

v frequency 

discrete frequency 

£ coordinate normal to shock wave 

p density 

p Q reference density, 1.292 kg/m3 

a Stefan -Boltzmann constant 

t time 

(f> angle between direction of propagation of radiation and y-axis 

\p azimuthal angle 

£2 solid angle 

Subscripts: 

b at body surface 

j,j at base point of bicharacteristic 

max maximum 


s 


immediately behind shock wave 



st 




at stagnation point 
w at body point 

OO in free stream 

Superscript: 

k number of time steps, t = k At and r = k At 

ANALYSIS 

The procedures which are used to solve the equations for time-dependent, inviscid 
flow at points within the shock layer, on the shock wave, and on the body surface are 
described herein. In addition, discussions of the thermodynamic and absorption coeffi- 
cient models for high-temperature air are presented, the method for evaluating the inte- 
gral terms in the energy equation is given, and the calculation of the starting solution is 
discussed. It should be noted that results for steady flow are obtained after many time 
steps when the time derivatives, which approach zero asymptotically, are sufficiently 
small. 


Calculations at Points Within Shock Layer 

Governing equations.- The natural coordinate system for formulating the problem of 
axisymmetric flow past a sphere is a spherical polar system with its origin located at the 
center of the sphere and its axis alined with the direction of flow in the free stream. 

When written in this coordinate system and in conservation form the governing equations 
are: 

Continuity: 


9(pr) 9(pur) 9(pv) 
9t + 9r + 90 


+ pu + pv cot 0 = 0 


Radial momentum: 

9 (pur) 8 [(p + pu 2 )r] 9(p UV ) 

+ — + 

9t 9r 90 




+ puv cot 0 = 0 


Tangential momentum: 


9(pvr) 9(puvr) 
9t 9r 


a(p + pv 2 ) 


90 


+ 2puv + pv 2 cot 0=0 


( 1 ) 

( 2 ) 


( 3 ) 
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Energy: 


a(pEr) a(puHr) 3(pvH) 
at + ar + qq 


+ puH + pvH cot 6 = rQ 


(4) 


In these equations, the radial and tangential velocity components, the density, and pres- 
sure are represented by u, v, p, and p, respectively. The quantities E and H are 
the total energy and enthalpy, respectively, and are given by 


E = e + i(u 2 + v 2 ) 

Ct 

H = h + i(u 2 + v 2 ) 

Ci 

where e and h are the internal energy and static enthalpy, respectively. The quan- 
tity Q in equation (4) is the net heat input to a volume element due to the emission and 
absorption of thermal energy and is discussed in a subsequent section. 

Equations (1) to (4) are indeterminate on the axis where the tangential component of 
velocity v vanishes and cot 0 is infinite. The indeterminate forms can be evaluated 
with the aid of THospitaTs rule to yield the following equations: 

Continuity: 


9(pr) + a(pur) | 2 3(py) 

9t 9r 90 


Radial momentum: 


a(pur) , 9 [(P + Pu 2 )r] | 
9t dr 


2 ^v) . (p _ p U 2) = o 

do 


(5) 


(6) 


Energy: 


8(pEr) + a(puH r) + 2 ^pyH) + = rQ (7) 

at ar de 

The tangential momentum equation is not needed at 9 = 0 because the tangential velocity 
component is zero on the axis. 

Equations (1) to (7) can be written in index notation as: 
aAi aBi dCi 

It +- a7 + (1 + M) ^F +Di + (1 " M)E i COt0 = F i (i = 1,2,3, 4) (8) 
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where 


and where 


fJL = 0 

M = 1 


( 0 / 0 ) 


( 0 = 0 ) 


( 9 ) 


= pr 

Ag = pur 

= pur 

Bg = (P + pu 2 )r 

o 

II 

Cg = PUV 

Dj = pu 

Dg = -P + p(U 2 - V 2 ) 

= pv 

Eg = pUV 

*1 

II 

o 

Fg=° 

Ag = pv r 

A 4 = pje + |(u 2 + v 2 )jr 

Bg = puvr 

B 4 = pujh + |(u 2 + v 2 )j : 

C 3 = p + pv 2 

C 4 = pv[h + |(u 2 + v 2 )] 

Dg = 2puv 

D 4 = pujh + |(u 2 + v 2 )J 

Eg = pv 2 

E 4 = pv jh + |(u 2 + v 2 )j 

o 

II 

CO 

Pu 

F 4 = Qr 


New independent variables are introduced as follows: 



0=0 


(Hb) 


Y = 


r - r b 
6 


(He) 


where 


6 = 6(0, t) = r s (0,t) - r b 


( 12 ) 


The quantity r g (0,t) is the radial coordinate of the shock wave for particular values of 
0 and t, and r b is the radius of the sphere. Equations (8) are written in terms of 
these independent variables as 



r aA i 

• w — - 

L 3Y 


+ (1 + M) 


a^actll 

30 3YJ f 


+ ( 1+ ^W 


+ D| + (1 - p.)E| cot 0 = 


(i = 1,2, 3, 4) (13) 


where the component of the shock velocity normal to the body is given by 


w = |f (14) 

The partial differential equations (13) are approximated with finite-difference equa- 
tions at points within the shock layer. These finite -difference equations and the known 
solution at time r are used to determine the density p, the internal energy e, and the 
velocity components u and v at time t + At. The pressure p must be determined 
from an equation of state as a function of p and e. The static enthalpy h is defined 
by the equation 

h = e + ^ (15) 

P 

In order to determine the equation of state, it is convenient to define the thermody- 
namic function y(p,e) as the ratio of static enthalpy to internal energy: 

y(p,e) = | (16) 

By use of equations (15) and (16), the required equation of state is given by 

p = pe(y - 1) (17) 
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It follows that 


e = -i_E 

y - 1 P 


h = — 2L-£ 
y - 1 P 




(18) 


Equations (17) and (18) are the same as the calorically perfect gas relations when y is 
replaced by the ratio of specific heats y. However, equations (17) and (18) are true for 
all gases, regardless of the thermodynamic state. 

Method of solution .- Let k, l, and m be the grid indices for the t, Y, and 
0 coordinates, respectively. The computational grid associated with these coordinates 
is shown in figure 1. As the figure indicates, a staggered grid is used. 



Figure 1.- Floating coordinate system. 
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The partial derivatives BAi/Bt in equations (13) are replaced by forward differ- 
ences as follows: 



1 

4 


(A-i) 


k 

l+l,m 


+ (Aj) 


k 

l-l,m 


+ (Aj) 


k 

l,m+l 



+ 0(A) 


(19) 


These difference expressions contain "artificial viscosity" error terms of first order 
which serve to stabilize the numerical calculations. The partial derivatives BB^/bY 
and 8Cj/90 at the point (k,l,m) are replaced by central differences of the form 



1 

2 AY 
1 

2 A0 


[®C.» 



+ Q(A 2 ) 



+ 0(A 2 ) 


(20a) 


(20b) 


The difference expressions for the derivatives BA^jBY and BC^jBY are similar to 
those for 9Bj/9Y. 

Since a staggered grid is used, the flow properties are not known at the point 
(k,l,m). Therefore, the undifferentiated terms and are replaced by the 
averages of the values of these quantities at the neighboring grid points: 


(Di) 


1, m 


, .k , k . ,k , k 

(Di) + (Di) + (Di) + (Di) 

1+1, m 1-1, m l,m+l l,m-l 


0(A 2 ) 


(21a) 



(Ei)f + (Ei)f + (Ei)f + (Ei)f 

1+1, m 1-1, m l,m+l l,m-l 


+ 0(A 2 ) 


(21b) 


The terms are all zero except F4 which is equal to rQ. The treatment of this 
term is discussed in the section "Radiative Transfer." The expressions in equations (19), 
(20), and (21) and the values of the terms Fi are substituted into equations (13), and the 
quantities 6, W, and 9 5/90 are evaluated at the point (k,l,m). All the terms in the 

k+ 1 

resulting finite -difference equations have a time index k except the quantities (Aj) 

l,m 

Therefore, values for A^ terms at time k + 1 can be evaluated if the solution at 
time k is known. 
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irj. 1 k+ 1 k+ T 

The velocity components u 7 and v 7 , the density p ? rn , and the internal 

energy e? H "* are determined from the quantities (Aj)^ + ^. Then the function is 

*> m £,m Z,m 

k+ 1 k+ 1 

evaluated in terms of p, m and e 7 The method by which this is done is discussed 

't, m Z,m J 


k+1 

later in the section "Thermodynamic Models." Values for the pressure p 7 and the 

C f XU 

k_L 1 

static enthalpy h^ m are determined from equations (17) and (18), respectively. 

It is seen from equations (10) that the radial coordinate r must be known at the 
point (k+1, l, m). The equation for r is 



+ Y6 


k+1 

m 


k+1 

Since the shock-layer thickness 5 m appears on the right-hand side of this equation, 
the solution at the shock at time k+1 must be determined first. 

The solutions at the boundary 0 = © max which lies downstream of both the sonic 
line and limiting characteristic are determined by extrapolating the flow properties lin- 
early. The errors incurred by the extrapolation do not affect the solution upstream of the 
boundary because the flow in the vicinity of the boundary is supersonic. 

Stability analysis .- A linear stability analysis of the finite -difference equations (13) 
has been performed by using the method of Von Neumann (ref. 13). The purpose of the 
resulting stability conditions is to insure that the magnitudes of the infinitesimal errors 
which are introduced at a given time step do not increase with time. Von Neumann showed 
that these errors are wavelike in nature. Although the quantities of interest are the sums 
of the errors affecting the various flow properties, the stability conditions are developed 
by considering the growth of only one error term for each flow property. 

The Von Neumann conditions for the finite -difference equations employed in this 
report are 


Ar < \fl 6 

AY _ 2 |u - YW| + a 


At <\^2 r b + Y6 
A0 “2 | v | + a 


where a is the speed of sound. The derivation of these inequalities is given in refer- 
ence 14. It should be noted that these conditions are more restrictive than the 
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Courant-Friedrichs-Lewy conditions (ref. 15 ), which do not include the factor \[ 2/2 on 
the right-hand side. 

Calculations at Points on the Shock Wave 

Governing equations .- An orthogonal curvilinear coordinate system 77, £ is estab- 
lished at the shock wave as shown in figure 2, where 77 and £ are the distances along 
and normal to the shock, respectively, 
at time t + At. The velocity com- 
ponents in the £ and 77 directions 
are U and V, respectively. The 
angle e between the normal to the 
shock and the radial line through the 
point being computed is determined 
from the relation 


tan e =^§ 
r b + 6 


( 22 ) 


The equations which must be 
solved at the shock wave are the 
Rankine-Hugoniot relations and one 
characteristic compatibility relation. 

If the subscripts °° and s repre- 
sent quantities in the free stream and 
immediately behind the shock, respec- 
tively, the Rankine-Hugoniot relations 
can be written as 



Figure 2 - Shock-oriented coordinate system. 


p s - s (” s -^) 2 = p »+ p ~(°~-^) 2 

l(u s . -5L-) 2 , ifu„ 

_ 1 P S 2 \ COS 6 / y^ _ 1 2 \ 


> ( 23 ) 


W 


v s = v co 


cos e , 


15 



The compatibility relation, which is a linear combination of the governing partial differ 
ential equations (ref. 16), is given by the equation 


pa^U - 1 (pav pa2 |V + v 2 e\ _ A. 
dr dT h^\ dr] 8rj dr]) h^ 


pa 2 U 


K + h? ? R Sin P) + pa2y N k C ° S 



+ 




(24) 


where the angle between the normal to the shock wave and the free-stream direction /3, 
the curvature of the shock wave K, and the scale factor for the 77 coordinate h^ are 
given by the equations 


/3(?7) = 0 - e 

K <”> ■ If 

h v (v,t) = 1 + K£ 

The perpendicular distance from the axis to the point 77, £ is given by the equation 


r rbv 

R07,£) = \ cos 73(77) d77 + £ sin 73(77) 

J 77=0 

where 77 = 0 at the axis. Two of the undifferentiated terms on the right-hand side of 
equation (24) are indeterminate on the axis where sin 73, V, and R are equal to zero. 
With the aid of 1' Hospital’s rule, the compatibility relation on the axis can be written as 


dr dr 


J-(pa 2 |X. + Kpa 2 U 

M dr l 


y 


- 1 + 


9 y 


9 log e e 


Q 


(25) 


It should be noted that the form of the energy equation which is used to obtain the 
compatibility relation is written as 


dp „2 d£ _ 

dT ‘ dr" 



9 y 

9 logp 


Q 


(26) 


where the speed of sound is given by the equation 


a 2 = e )(y - 1) 




(27) 




dy 

8 log p 



Equations (26) and (27) are derived in appendix A. 

The compatibility relation is integrated along the line of intersection of the 
£,t plane which passes through the shock point at time r + At and the characteristic 
surface which passes through the shock point and intersects the f], 4 plane behind the 
shock wave at time r. This line of intersection has the slope 


= U + a (28) 

dT 


in the £,t plane. 

Method of solution .- The solution is known completely at time t and at all pre- 
vious times. A first approximation to the shock-wave geometry at time r + At is 
made by using the relation 

6(0, t + At) = 6(0, t - At) + 2W(0,t - At) At 

and the component of the shock velocity normal to the body is approximated by 

W(0,t + At) = W(0,t - At) 

The derivative 36/80 is computed from the equation 


85 _ 5(0 + 2 A0, t + At) - 5(0 - 2 A0,t + At) /ggs 

80 4 A0 1 ^ 

and the angle e is determined from equation (22). Tentative values for the velocity 
components U s and V s and the equilibrium thermodynamic properties p g , p g , e s , 
h s , a s , and y g at time t + At are obtained by using the approximate value for 
W(0 ,t + At) and iterating the Rankine-Hugoniot relations (23), equations (18) and (27), 
and the equation for y(p,e). 

At this point, the location of the line along which the compatibility relation is to be 
integrated is determined. The slope of this line, which lies in the £,t plane and passes 
through the point s, is given by equation (28). The £,t plane is shown in figure 3. 
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Figure 3.- Location of line along which compatibility equation is integrated. 

The first estimate of the location of the point j where the line intersects the rj, £ plane 
at time r is obtained by using the relation 


= -(U + a) g At 

where £ is equal to zero at the shock, and its positive direction is outward into the free 
stream, as shown in figure 3. The mesh spacing At is always positive. The values of 
the flow properties at point j are determined by interpolation, and is computed 
again from the relation 

= -|[(U+a) s + (U + a)i]AT 

so that a new location for point j is determined. This process is repeated until the 
successive values for converge. 

The compatibility relation given by equation (24) or equation (25) is integrated 
from point j to point s. In order to perform this integration, the right-hand side of 
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equation (24) or (25) and the term pa on the left-hand side are evaluated at s and j, 
and the results are averaged. The integrated compatibility relation, the Rankine-Hugoniot 
relations, and the thermodynamic relations are solved iteratively for the velocity com- 
ponents, the thermodynamic properties, and the shock velocity component W at the 
shock point at time r + At. 

The next approximation of the shock-wave geometry at time r + At is determined 
by using the relation 

6(0, t + At) = 6(0, t - At) + [w(0,t + At) + W(0,t - At)] At (30) 

The process given by equations (29) to (30) is repeated until the successive solutions 
converge. 


Calculations at Points on Body Surface 

Governing equations .- The transient spherical polar coordinate system r,0,t is 
used to obtain solutions at the surface of the sphere. The differential equations which are 
solved are the energy equation (26), the tangential momentum equation, and a character- 
istic compatibility relation. Both the energy equation and the momentum equation are 
integrated along the body streamline, which has the slope 


d0 _ v_ 
dt r b 


(31) 


in the 6 , t surface. Equation (26) is used in the form given, whereas the momentum 
equation is used in the form 


dv _ 1 8p 

dt pr 9 0 


(32) 


The compatibility relation is integrated along the line of intersection of the r,t plane 
(0 = Constant) which passes through the body point where the flow is being computed at 
time t + At and the characteristic surface which passes through the body point and 
intersects the r, 0 plane within the shock layer at time t. This line of intersection has 
the slope 


— = u - a (33) 

dt 

in the r,t plane. The characteristic compatibility relation to be used at the body sur- 
face is 
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( 34 ) 


_ pa — = 
dt p dt 


= i/pav |H - 

rr 9 0 


pa * F? " v ig) " p( 2 P a2u + P a2y cot 9 + pav 2 ) 


, 2 9v _ v 9£) 


1 + 


3y 


i 8 lQ g^ e > 


Q 


for points off the axis of symmetry and 


<>E - pa |U - 2 ( pa 2 2 

dt dt rl be 


(pa 2 §I + pa2 u ) + p- 1+ (_g_) 


Q 


(35) 


at the stagnation point. 

Method of solution .- A first approximation to the solution at time t + At (which 
is the same as t + At since t = t) is obtained by integrating equations (13), the partial 
differential equations which are used at points in the flow field. Let k, 2 = 1, and m 
be the indices of the body point w for the t, Y, and 0 coordinates, respectively. 
The partial derivatives 9A^/9 t and 9Bj/9Y in equations (13) are approximated by 
finite-difference expressions of the form 
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and 
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AY 



k 
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+ O(A) 


respectively. Since Y = 0 at the body surface, the partial derivatives 9Aj/9Y and 
9Cj/9Y need not be evaluated. The derivatives 9Ci/90 are determined from equa- 
tion (20b); and the quantities and are evaluated at the points k, 1 = 1, 
m + 1, and k, 1 = 1, m - 1, and the results are averaged. 

The location of the line along which the compatibility relation is integrated is deter- 
mined from equation (33). The first estimate of the location of the point j where this 
line intersects the r, 6 plane at time t is obtained from the relation 


r j = r b + At 


(36) 
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The values of the flow properties at point j are determined by interpolation, and rj 
is computed again from the equation 

r j = r b + |[ a w + ( a - u )j] At 

This process is repeated until the successive values for rj converge. The compatibility 
relation given by equation (34) or equation (35) is integrated from point j to point w 
in the same manner as described for the shock wave. Because u = 0 at the surface, 
the integration of the compatibility relation yields the new value for the pressure p w at 
time t + At. 

The location of the line along which the energy equation (26) and the tangential 
momentum equation (32) are integrated is determined from equation (31). When the time 
index k for the time plane t + At is even, this line extends from the body point w 
at time t + At to a body point 3 at time t - At; when k is odd, the line extends from 
point w to a body point 3 at time t. In other words, point 3 is always located on a 
time plane with an even index k. It can be seen from figure 1 that the set of body points 
at which solutions are obtained when k is even includes the stagnation point. This pro- 
cedure of integrating equations (26) and (32) over one and two time intervals alternately 
has been found to yield better results than those which were obtained by always integrating 
over one time interval. The first estimate of the location of point 3 is determined from 
the relation 


or 


e V 
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depending upon whether k is even or odd, respectively. Successive estimates of 9 ^ 
are determined from the relation 
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The integration of equations (26) and (32) from point j to point w is accomplished 
by evaluating the right-hand sides of these equations at these points and averaging the 
results. The term a 2 on the left-hand side of equation (26) is evaluated in the same 
manner. Equation (26) yields the new value of the density p w since the new value for 
p w is known, and equation (32) yields the tangential velocity v w . 

The process described, starting with equation (36), is repeated until the successive 
solutions converge. 


Thermodynamic Models 

Two thermodynamic models are used in this report: an equilibrium air model and 
a perfect gas model. Both models use equations (17) and (18) to determine the pres- 
sure p and enthalpy h as functions of the density p and internal energy e. The 
ratio y = h/e is assumed to be constant for the perfect gas model, whereas curve fits 
are used to evaluate it for the equilibrium air model. 



(a) Ratio of enthalpy to internal energy. Data from references 17 and 18. 
Figure 4.- Curve fits of equilibrium air thermodynamic data. 
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Equilibrium air model .- A curve fit is used to describe the variation of the ratio 
y = h/e as a function of p/p Q and e/RT 0 as shown in figure 4(a). The data used to 
construct this fit were obtained from the results of Allison (ref. 17) and Browne (ref. 18). 
It can be seen from the figure that for a given density, the variation of y with the loga- 
rithm of e/RT 0 can be characterized by straight-line segments over most of the range 
1 S e/RT 0 § 1500. Transitions from one straight-line segment to another occur for values 
of e/RT 0 near 6, 25, 70, and 400. The figure shows that the transition in the vicinity 
of e/RT 0 = 25 occurs rapidly, and it is assumed in this report that this transition can 
be approximated by an abrupt change from one straight-line segment to another. The 
same assumption is made with respect to the transition in the vicinity of e/RT 0 = 6, 
although this transition does not occur as quickly. Transitions in the vicinity of 
e/RT 0 = 70 and 400 are smooth, and they are approximated with odd transition functions 
(ref. 19). The curve fit for y which is used in this report is given by equations (Bl) 
to (B4). 

The curve fit which is used to determine the temperature in terms of the pressure 
and density is compared with the data of Browne in figure 4(b). This property is needed 



(b) Thermal equation of state. Data from reference 18. 
Figure 4.- Continued. 


23 



only in connection with the computation of the heat flux and the divergence of the heat flux 
for the radiation problem, which are discussed in a subsequent section. The ranges of 
interest for the calculations in this report are 


10-3 s£_< 10" 1 
^o 

8 000° K ^ T ^ 16 000° K 

The data of Allison and Browne show that for sufficiently large values of the temperature 
and small values of the density within this range, the temperature can be predicted very 
accurately with a power -law fit of the form T « p-432/ p-392. Similarly, for sufficiently 
small values of the temperature and large values of the density, the temperature can be 
predicted with a power-law fit of the form T p .810/^.794. j n report, an even tran- 
sition function (ref. 19) is used to effect the change from one power -law fit to the other. 
The expressions which are used in this report to determine the temperature in terms of 
the pressure and density are given by equations (B5) and (B6). 

The speed of sound is used in the determination of the flow properties at the shock 
wave and the body surface. Equation (27) is used to determine the speed of sound in 
terms of p and e in this report. The derivatives of y which appear in equation (27) 
are obtained by differentiating the expressions for y in appendix B. The results of this 
procedure are compared with the data of Browne in figure 4(c). 

Perfect gas model .- This model is formulated subject to the assumptions that the 
ratio y - h/e is constant and that the temperature is given by the equation 


T = 




(37) 


where the exponents R and S differ from one. It is shown in appendix C that if the 
entropy S is to be a conserved quantity, the ratio y and the exponents R and S 
must be related by the equation 


1 - S = y(l - R) 


(38) 


In other words, only two of the three parameters y, R, and S can be chosen indepen- 
dently. It is shown also in appendix C that the speed of sound for this model is given by 
the relation 

a 2 = 

P 
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e/RT 

o 


(c) Speed of sound. Data from reference 18. 
Figure 4.- Concluded. 


Radiative Transfer 

In order to determine the amount of thermal energy which is absorbed by a fluid 
element, it is assumed that the radiating shock layer can be approximated by an infinite 
slab with the thickness of the shock layer at the point under consideration. (See fig. 5.) 
The distributions of the thermodynamic properties across the slab are the same as those 
along a normal to the surface through the point. 

Both frequency-independent (gray) and frequency -dependent (nongray) absorption 
coefficient models are treated. The Planck mean absorption coefficient is used for the 
gray model, whereas a two-step coefficient is used for the nongray model. One step lies 
in the vacuum ultraviolet portion of the spectrum (400A IAS 1130A) and the other covers 
the near ultraviolet, visible, and infrared portions of the spectrum (1130A = A < «>). A 
Planck mean averaging procedure is used to determine the absorption coefficients for the 
two spectral bands. 

One -dimensional radiating shock layer model .- The equation for the net addition 
of heat to a volume element due to monochromatic radiation can be obtained from 
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Slab 



Figure 5.- Schematic of siab used for calculation of heat-flux 
and net-heaHnput terms. 


equation (D14). This equation is an expression for the net monochromatic heat input per 
unit volume for a one -dimensional slab. It was derived subject to the assumptions that 
the body surface neither reflects nor emits radiation, and that the gas in the free stream 
neither emits nor absorbs. The quantities I^(y^) and Ij,(y s ) therefore have values of 
zero. When equation (D14) is substituted into equation (D15), the following expression 
for Q(y) is obtained: 


Q(y)=-2 a I/ ( 5 r)^B l ,(jr b )E 2 [A J/ (y b> y)] + B y (y s )E 2 [A„(y,y s )] 


r»s=y i dB„ 

+ 7— r — ^ d ^ E< 
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where Bj, is the Planck function 
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E n (z) is an exponential integral 
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(39) 


p(. o=°° -o>z 
E n (z) = \ 
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(40) 


and Aj,(r,s) is the optical depth 

A„(r,s) = f oi u (t)dt 

J t=r 

Consider an absorption coefficient model with M steps: 

<*v = ^ (^N = v = "n+1 ;N = 1 * 2 >* • *’ M ) ( 41 ) 


The integral of the Planck function over a frequency interval is given by the equation 
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where 
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The net heat input for an M step nongray absorption coefficient model is 
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(44) 


Equation (44) is used to determine the quantity Q in equations (4) and (7). 
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Equations (D7) and (D8) can be combined to obtain an equation for the heat flux at 
the body surface q(yb): 


q(yt>) = - 2 * 




The expression for the heat flux at the surface for an infinite slab and an M step absorp- 
tion coefficient model is written as 


q(Yb) = ~ 2 ° 



T 4 (s) 




(45) 


Equations (44) and (45) are integrated numerically. The absorption coefficients 
are used in tabular form as functions of the density p and the temperature T. 

Absorption coefficient m odels.- The one- and two-step models which are used in 
this report are constructed from the more complex models of Olstad. (See ref. 5.) The 
sources of radiation which are treated in that reference include the following groups: 

(1) atomic lines, (2) N + and 0 + electron recombination in the vacuum ultraviolet portion 
of the spectrum, and (3) molecular band systems, free-free electron transitions, the 
photoionization of 0, and electron recombination in the visible and infrared portions of 
the spectrum. 

Olstad determines the absorption coefficients for electron recombination radiation 
in the vacuum ultraviolet portion of the spectrum from the absorption cross sections of 
Hahne (ref. 20). The data of Allen (refs. 21 and 22) are used to compute the absorption 
coefficients for the radiation sources listed in groups 1 and 3. 

The step of the nongray model for the vacuum ultraviolet portion of the spectrum 
(400A ^ = 1130A) accounts for the electron recombination radiation in that spectral 
region. The radiation from this part of the spectrum is optically thick for all the shock 
layers treated in this report except those associated with spheres with radii of 1 foot or 
less. 

The second step of the nongray model accounts for the radiation due to the processes 
listed in group 3, the atomic lines above 2000A, and the wings of the lines below 2000A. 

The radiation due to these processes occurs in the near ultraviolet, visible, and infrared 
portions of the spectrum (1130A = X ^ °°) and is optically thin for the cases treated in this 
report. 
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(a) Gray model. 

Figure 6.- Absorption coefficients for equilibrium air. 


The present model does not account for the radiation from the centers of the lines 
below 2000 A, which is optically thick. In order to take it into account, a model with more 
than two steps would have to be employed. The errors encountered by neglecting the 
radiation from the line centers are probably not severe since this radiation tends to be 
absorbed quickly. 

The gray absorption coefficient model which is used in this report is shown graphi- 
cally in figure 6(a), and the nongray model is shown in figure 6(b). 

Starting Solution 

The starting solution is constructed by specifying the shock-wave shape and surface 
pressure distribution. Since the shock wave is assumed to be stationary initially, the flow 
properties at points on the shock wave can be determined from the Rankine-Hugoniot rela- 
tions. The properties at points on the body surface are determined by using the specified 
pressure distribution and assuming that the entropy and total enthalpy are conserved along 
the body streamline. The pressure and the velocity components are determined at points 


29 




8 000 10 000 12 000 14 000 16 000 

T , °K 

(b) Nongray model. 

Figure 6.- Concluded. 

within the shock layer by interpolating the values at the shock wave and body surface 
linearly along lines normal to the body surface. The density and enthalpy profiles along 
the stagnation streamline are determined with a solution similar to that used by Goulard 
in reference 1. This solution is developed in appendix E. Similar profiles are used 
along the other lines normal to the body surface. The values of the stagnation -point 
enthalpy which are used for the starting solution are obtained by linearly extrapolating 
the values at the two neighboring points on the stagnation streamline. (See appendix E.) 
The correct values for the cases treated in this report are much lower than the extrapo- 
lated values since a radiation -cooled layer exists at the stagnation point which is much 
thinner than any practical value of the mesh spacings which can be used. The use of the 
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correct value of the stagnation -point enthalpy would lead to an erroneous broadening of 
the radiation-cooled layer over several mesh spacings since the first-order computa- 
tional method which is used in this report tends to reduce steep gradients in the flow 
property profiles such as those associated with the radiation -cooled layer. It has been 
found that better solutions for the flow field are obtained when the extrapolated value of 
the stagnation-point enthalpy is used, although the radiation-cooled layer is eliminated by 
this process. This procedure does not affect the value for the stagnation -point heat flux 
appreciably since the radiation -cooled layer is so thin. It should be noted that the pres- 
ent computational method adjusts the flow properties at the stagnation point much more 
slowly than elsewhere. Therefore, the starting solution for the stagnation streamline, 
and hence the extrapolated starting solution for the stagnation point, must be determined 
as accurately as possible. 


RESULTS 

In this section, calculations which are determined with the present method for both 
adiabatic and nonadiabatic flow are compared with those obtained with other techniques. 
The stagnation streamline enthalpy profiles and the shock-layer thickness distributions 
for spheres with different radii and the same flight conditions are shown, and the gray 
and the nongray surface heat-flux distributions for these spheres are compared. In addi- 
tion, the effects of sphere size, free-stream velocity, and altitude on the nongray surface 
heating distributions are investigated. The free-stream thermodynamic properties which 
are used in this report are based on the atmospheric model of reference 23. 

Illustrative Calculation for Adiabatic Flow 

In order to demonstrate the accuracy of the present computational method and 
thermodynamic model, an illustrative calculation for adiabatic flow is made. The com- 
putation is made for the flow of equilibrium air past a sphere traveling at 45 000 ft/sec 
(13.73 km/sec) at an altitude of 200 000 ft (60.96 km) in the atmosphere of the earth. 

The mesh spacings which are used are AY =1/6 and A0 = 3.875°. 

The results of the present method are compared in figure 7 with those of Lomax 
and Inouye (ref. 24), who use the inverse method. It is seen in figure 7(a) that the results 
of the two methods for the shock-wave and sonic -line locations are in close agreement. 
The present method yields a value of 6/rjj = 0.0467 for the stagnation streamline shock 
detachment distance, whereas the inverse method gives the value 6/r^ = 0.0479. The 
results for the tangential velocity and density distributions at the surface are shown in 
figure 7(b). The distributions of these flow properties along a line normal to the body 
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(a) Shock and sonic line locations. 

Figure 7.- Adiabatic flow about a sphere traveling at 45 000 ft/ sec (13.73 km/sec) 
at an altitude of 200 000 ft (60.96 km). 

surface at 0 = 38.75° are given in figure 7(c). It should be noted that the present 
method computes the density less accurately than any of the other flow properties. 

Sample Calculations for Nonadiabatic Flow 

The gray radiation model is used for these sample calculations. The method used 
to obtain the initial solutions has been discussed previously. Equations (E27) and (E29) 
are used to determine the initial stagnation streamline enthalpy distribution. A one-step 
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absorption coefficient which covers the spectral range 400A = X < 00 is used for 
the gray radiation model. The power -law expression, which is obtained by curve 
fitting the data presented in figure 6(a) over the ranges 10“3 s < 10“ 1 and 
10 000° K ^ T ^ 16 000° K, is given by the equation ° 


/o \1.179 / T v 
Q!p= 0.3060 £-) -Lr) 

P \ p o Uo 4 / 


.6.964 


cm 


-1 


The dependence of temperature on enthalpy is assumed to be 

T cc h- 392 


The constant C in equation (E27) is set equal to zero since a one-step absorption coef- 
ficient is being considered. 

Effect of mesh spacin g size.- Sample calculations are made for a sphere with a 
radius of 5 feet (1.524 m) traveling at a speed of 50 000 ft/sec (15.24 km/sec) at an alti- 
tude of 190 000 ft (57.91 km) in the atmosphere of the earth. 

The influence of the mesh spacing size on the enthalpy distribution is shown in fig- 
ure 8. Calculations are presented in figure 8(a) for AY = 1/14 and A 0 = 1°, 2°, 4°, 
and 6°; the solutions for the enthalpy distribution and the shock -detachment distance are 
seen to converge as A0 is decreased. In figure 8(b) results are shown for A© = 2° 
and AY = 1/6, 1/10, 1/14, and 1/18; the results are seen to converge as AY is 
decreased, although the convergence is not as strong for AY as for A0. 

The effects of the size of the mesh spacings A0 and AY on the surface heat- 
flux distribution are shown in figures 9(a) and 9(b), respectively. There is a strong ten- 
dency for these distributions to converge as the size of the mesh spacings is reduced. 

The results for the distributions are in close agreement for AY = 1/14 and A© = 1° 
and 2°. For A© = 2° and AY = 1/14 and 1/18, the results coincide for 0 S 30° and 
show close agreement for © > 30°. 

On the basis of the results presented in figures 8 and 9, the mesh spacings for the 
subsequent calculations are chosen to be AY = 1/14 and A© = 1°. 

Comparisons with other methods .- The present method has been used to calculate 
two cases which are presented by Olstad in reference 4. Olstad employs the optically 
thin approximation and uses the Poincare-Lighthill-Kuo method to make the calculations 
for these cases. 

The perfect-gas thermodynamic model is used to calculate the results of the pres- 
ent method. The free-stream pressure is set equal to zero and the quantity y = h/e 
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is chosen so that the density at the shock wave is — = 0.06, the value used in refer - 

Ps 

ence 4. The dependence of temperature on the density and pressure is of the form 

d .601 

T oc c 

p.550 

The density exponent S = 0.55 is the same as the enthalpy exponent used in reference 4, 
and the pressure exponent R was computed from equation (38). The dependence of the 
absorption coefficient on the density and temperature is taken to be 


a-p ~ p 


.661 t 6.656 


Therefore, the absorption coefficient varies as h^ when the pressure is constant. This 
variation was used in reference 4 for the cases which are being treated. 


Let the energy depletion parameter T be defined by the equation 


T = 


4a P,s^A gT s 

P coVoohg 


4 


(46) 


where 6^ is the shock-detachment distance for adiabatic flow. The gases which are 
treated have the parameters 0 !p g 6^ = 0.3, T = 0.02 and Op g 6^ =1, T = 2. 


The results of the present method for the stagnation streamline enthalpy profiles 
are compared with those of Olstad in figure 10. The agreement is reasonably good for 
both cases. The results of the present method for the stagnation -point heat fluxes are 
q„t 

— _ ■ = 0.0032 and 0.0632 whereas the values presented in reference 4 are 0.0028 and 

p V 3 

0.0574. The present method gives values of 6 = 0.9946^ and 0.7976^ for the shock- 
detachment distances for the two cases, whereas reference 4 gives values of 0.9826^ 
and 0.7326» . 


The results of the present method for the surface heat -flux distribution are com- 
pared with those of Wilson and Hoshizaki (ref. 7) in figure 11. The case which is treated 
is a sphere with a radius of 5 feet (1.524 m) traveling at 50 000 ft/sec (15.24 km/sec) at 
an altitude of 190 000 ft (57.91 km) in the atmosphere of the earth. The gray radiation 
model is used. Two distributions from reference 7 are shown in the figure because the 
authors use upper bound and lower bound absorption coefficient models in that reference. 
It is seen that the heat -flux distribution computed by the present method is very similar 
to the distribution of reference 7. 
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Calculations for Spheres With Different Radii and 
Different Flight Conditions 

The thickness of the shock layer is roughly proportional to the radius of the body. 
Therefore, the body size governs the optical depth of the shock layer. The flight velocity 
determines the amount of energy available to the flow field, whereas the velocity and the 
free-stream conditions together determine the magnitudes of the flow properties behind 
the shock wave. 

Altitude, 190 000 feet (57.91 km); velocity, 50 000 ft/ sec (15.24 km/sec) .- Results 
for spheres traveling in the atmosphere of the earth at an altitude of 190 000 feet 
(57.91 km) at a velocity of 50 000 ft/sec (15.24 km/sec) are presented in this section. 

The effects of body size on the stagnation streamline enthalpy profile, the shock -layer 
thickness, and the surface heat -flux distributions are investigated. Spheres with the 
radii of 3 inches (0.0762 meter), 1 foot (0.3048 meter), 5 feet (1.524 meters), and 
15 feet (4.572 meters) are treated. 

The initial solution for the gray cases has been discussed previously. A two-step 
absorption coefficient model is used in connection with equations (E27) and (E29) to deter- 
mine the initial solution for the nongray cases. Again, the temperature is assumed to 
vary as h-392_ The data presented in figure 6(b) are used to obtain the power-law 
expressions for the absorption coefficient model. The coefficient for the step which 
covers the spectral range 1130A g X < °° is written as 



where the temperature T is given in °K. For the cases r^ = 3 inches and r^ = 1 foot, 
the second step covers the spectral range 400A = X = 1130A. The power-law absorption 
coefficient for this range is 



For the case with r^ = 5 feet, the second step covers the spectral range 
1020A = X = 1130A. For this case, the radiation for A < 1020A is optically thick. It is 
assumed that for these conditions this radiation, which is absorbed by the flow field very 
close to the point at which it is emitted, does not affect the structure of the flow field 
appreciably and can be neglected. The power -law absorption coefficient for the range 
1020A = A = 1130A is given by the expression 
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For the case with = 15 feet, the radiation for A < 1130A is optically thick and is 
neglected. Therefore, only one step is used for this case. 



Figure 12.- Effect of size on gray and nongray stagnation-streamline enthalpy distri- 
butions for spheres traveling at 50 000 ft/ sec (15.24 km/sec) at an altitude of 
190 000 ft (57.91 km). 

In figure 12, the enthalpy distributions along the stagnation streamline are shown 
for the various cases. It is seen that the effects of radiation on the enthalpy are more 
pronounced for the gray model than for the nongray model. As the radius of the body is 
increased, the differences in the gray and nongray distributions increase. 

The variation of the shock-layer thickness with the angle 6 is shown in figure 13 
for the different cases. One effect of body size is an increase in the deviation of the 
shock-layer thickness for the nonadiabatic cases from that for the adiabatic case. When 
r^ = 3 inches, there is virtually no difference in the shock-layer thickness for the gray 
and the nongray cases. The nongray shock layers are thicker than the gray ones for the 
cases with r^ = 1 foot and r^ = 15 feet. For the case with = 5 feet, the gray shock 
layer has virtually the same thickness as the nongray one in the stagnation region, and it 
is thicker than the nongray one for larger values of the angle 9. 
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Figure 13.- Effect of size on gray and nongray shock-layer thickness distributions 
for spheres traveling at 50 000 ft/ sec (15.24 km/sec) at an altitude of 190 000 ft 
(57.91 km). 


The shock -detachment distance is governed largely by the rate at which mass is 
transferred out of the stagnation region. Along the stagnation streamline, this rate is 
proportional to the product of the density and the tangential velocity gradient. The den- 
sity is greater for the gray radiation model than for the nongray model at a given point 
in the stagnation region because the density is inversely proportional to the enthalpy and 
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the enthalpy for the gray model is less than that for the nongray model, as is shown in 
figure 12. For the case with r^, = 5 feet, the tangential velocity gradient for the nongray 
model is greater than that for the gray model by an amount which is just enough to com- 
pensate for the difference in the density. Therefore, the shock-detachment distances for 
the two radiation models are the same for r^ = 5 feet. 

In figure 14, the surface heat -flux distributions for the gray and the nongray radia- 
tion models are compared. Values for the stagnation -point heat fluxes are given subse- 
quently. It is seen that differences in the distributions for the two models are small. 
There is a slight tendency for the differences to increase as the body radius is increased. 

The gray and the nongray nondimensional surface heat -flux distributions presented 
in figure 14 are compared in figures 15(a) and 15(b), respectively, to determine the effect 
of body size on the distribution shapes. It is seen that this effect is slight for both radia- 
tion models. This result is particularly surprising for the nongray model since the rela- 
tive contributions from the different parts of the spectrum vary considerably as the body 
radius varies. (See fig. 16.) For the case with r^ = 3 inches, the heat flux comes pre- 
dominantly from that part of the spectrum with A < 1130A. The contributions of the two 
parts of the spectrum are approximately the same for the case with r^ = 1 foot. The 
heat flux from that portion of the spectrum with A < 1130A is only about 10 percent of 
the total for the case with r^ = 5 feet. For the case with r^ = 15 feet, the ultraviolet 
portion of the spectrum contributes less than 1 percent of the total heat flux. Therefore, 
a spectral breakdown of this case is not shown in figure 16. 

Altitude, 140 00 0 feet (42.67 km), ^velocity, 34 000 ft/sec (10.36 km/sec).- Results 
for spheres traveling in the atmosphere of the earth at an altitude of 140 000 feet 
(42.67 km) at a velocity of 34 000 ft/sec (10.36 km/sec) are presented in this section. 
For the flight conditions used in the previous section, the temperature in the hotter por- 
tions of the shock layers varies with pressure and density as p .432^.392. hotter 

portion of the shock layer provides most of the radiant heat which strikes the body sur- 
face. Since the heat flux depends strongly on the temperature, it might be expected that 
different flight conditions would yield different heat -flux distributions. The flight condi- 
tions which are used in this section are chosen so that the temperature in the shock layer 
varies as p .810^.794. The nongray surface heat flux is calculated for spheres with 
radii of 3 inches and 5 feet. 

In order to determine the initial solution with equations (E27) and (E29), it is 
assumed that the temperature varies as h*^4 # Equation (47) is used to obtain the 
absorption coefficient for the spectral range 11 30 A = A < °° for the cases with 
rjj = 3 inches and 5 feet. For the case with r^ = 3 inches, it is assumed that the second 
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(c) r b = 5 ft (1.524 m). 

Figure 16.- Concluded. 

step covers the spectral range 852 A = X ^ 1130A. The absorption coefficient for this 
range is 


«2 


68.15/ -2- 

[Pr, 


1.182 


&) 


.329 


cm 


-1 


The radiation for X < 852A is optically thick and is neglected. For the case with 
rjj = 5 feet, the radiation for X < 1130A is neglected. 

The results for the total nondimensional nongray surface heat-flux distributions for 
the two cases are shown in figure 17. For the present set of flight conditions, as with the 
previous set, the distributions are insensitive to body sizes. However, as with the pre- 
vious flight conditions, the contributions to the heat flux from the portions of the spectrum 
with X < 1130A and X > 1130A differ for the two cases. A spectral breakdown for the 
case with r^ = 3 inches is shown in figure 18. A breakdown for the case with r^ = 5 feet 
is not given because more than 99 percent of the heat comes from that portion of the spec- 
trum with X > 1130A. 


Effect of flight conditions on surface heat -flux distribution.- The nongray surface 
heat -flux distributions for the two different sets of flight conditions for spheres with radii 
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6, deg 


Figure 19.- Approximate analytic solution for nongray surface heat-flux distribution 
for a sphere with a radius of 5 ft (1.524 m). 


of 5 feet are compared in figure 19. Since the variations of the pressure, the density, 
and the shock-layer thickness with the angle 0 are very nearly the same for the two 
cases, the difference in the heat-flux distributions must be due to the differences in the 
dependence of the temperature on the pressure and density and in the dependence of the 
absorption coefficient on the density and temperature. In order to demonstrate these 
effects, an approximate analytic solution for the surface heat -flux distribution is devel- 
oped in appendix F. In figure 19 the approximate results for the nondimensional heating 
distributions which are determined from equation (F9) show the same trends as the 
numerical results. The approximate solution accounts only for the radiation from the 
optically thin portion of the spectrum with A > 1130A since the radiation for A < 1130A 
is optically thick and is absorbed in the flow field. It should be mentioned that the approx- 
imate results for the stagnation -point heat flux which are determined from equation (F10) 
are greater than those which are determined numerically. This condition occurs because 
the analytic solution accounts for absorption only in an approximate fashion. For the case 
with V m = 50 000 ft/sec (15.24 km/sec) and h = 190 000 ft (57.91 km), equation (F10) 
yields a value of 0.660 for the nondimensional stagnation -point heat flux q s ^2a!^ g 6(0)crT s ^ 
for A > 1130A. The numerical value for the heat flux for this case is 0.527. For the 
case with = 34 000 ft/sec (10.36 km/sec) and ii = 140 000 ft (42.67 km), the analytic 
and numerical results for the nondimensional stagnation -point heat flux for A > 1130A 
are 0.432 and 0.353, respectively. 
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The exponents a-^ and d in equation (F9) are assigned the values 1.2 and 1.3, 
respectively, for both cases. For the case with = 50 000 ft/sec (15.24 km/sec) 
and ii = 190 000 ft (57.91 km), the exponents bj, R, and S and the ratios y, 

6(0)/^, and P s t/Ps have the values 2.3, 0.432, 0.392, 1.129, 0.0404, and 1.87, respec- 
tively. The respective values for these parameters for the case with V M = 34 000 ft/sec 
(10.36 km/sec) and h = 140 000 ft (42.67 km) are 6.9, 0.810, 0.794, 1.147, 0.051, 
and 1.29. 

Both the altitude and the velocity differ considerably for the two cases which are 
compared in figure 19. In order to determine the individual effects of altitude and veloc- 
ity, additional calculations were made for the nongray radiation model and a body with a 
radius of 5 feet. The results for the total- surface heat-flux distributions for these cal- 
culations and those presented in figure 19 are compared in figure 20. For a given flight 
velocity, the effect of altitude on the surface heat-flux distribution is small. The velocity 
influences the heat -flux distribution chiefly through the equation of state, which governs 
the dependence of the temperature on the pressure and density. As stated previously, 
this dependence is p. 432^.392 and p.810/p.794 } respectively, for velocities of 
50 000 ft/sec (15.24 km/sec) and 34 000 ft/sec (10.36 km/sec). For velocities of 
41 000 ft/sec (12.50 km/sec), the dependence of the temperature on the pressure and 
density varies from approximately p-432/ p-392 on the axis to p .810^.794 f or large 
values of the angle 0. 



Figure 20.- Effect of flight speed and altitude on nongray surface heat-flux distribution 
for a sphere with a radius of 5 ft (1.524 m). 
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Although the shapes of the total nongray heat -flux distributions are very similar 
for a given flight velocity and different altitudes, the relative contributions from the two 
portions spectrum to the heat flux vary with altitude. For a velocity of 50 000 ft/sec 
(15.24 km/sec), the contributions to the heat flux from the optically thin portion of the 
spectrum with X > 1130A are approximately 90 percent and 60 percent, respectively, 
for altitudes of 190 000 ft (57.91 km) and 220 000 ft (67.06 km). The optically thin por- 
tion of the spectrum contributes approximately 100 percent and 65 percent of the surface 
radiative heat flux for a flight velocity of 41 000 ft/sec and altitudes of 140 000 ft and 
220 000 ft. The respective contributions for a velocity of 34 000 ft/sec and altitudes of 
140 000 ft and 190 000 ft are approximately 100 percent and 85 percent. 

The effect of the shape of the nondimensionalized heat-flux distribution on the total 
radiative heating to the body can be seen by comparing the ratio of the integral of the 
heat flux over the surface area with 0 2 44° to the product of this area and the 
stagnation -point heat flux for the cases shown in figure 20. This ratio, which is written 
as 


p0=44° 

'0=0 


\ q(0)sin0 d© 

Jo. - 


p0=44° 

q et \ sin 0 d0 

Sl J 0=O 


has values of 0.32, 0.48, and 0.63, respectively, for flight velocities of 34 000 ft/sec, 

41 000 ft/sec, and 50 000 ft/sec. These numbers demonstrate that comparisons of 
stagnation -point heating rates may not yield correct comparisons of the total heating 
rates for an entry vehicle. This observation has been made for the case of transparent 
shock layers by several writers. (For example, see ref. 25.) 

Stagnation -point heat f lux.- The present results for the gray stagnation -point 
heating rates are compared with those of Thomas (ref. 2) and Olstad (ref. 4) in figure 21. 

4 

These rates are made nondimensional by dividing by the factor 2a-p s 6 A oT s which is 

the gray heat flux obtained from an isothermal slab with the thickness of the adiabatic 
shock layer when absorption is neglected. The quantity against which the rates are 
plotted is the gray energy depletion parameter r, which is defined by equation (46). The 
results of Thomas and Olstad have been interpolated for a flight velocity of 50 000 ft/sec 
(15.24 km/sec) and an altitude of 190 000 ft (57.91 km) in the atmosphere of the earth. 

It is seen that the present results agree more nearly with those of reference 2 than with 
those of reference 4. 

The results fcr the nongray stagnation-point heating rates are given in table I. 
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Figure 21.- Gray stagnation-point heat-flux distribution for spheres traveling at 50 000 ft/ sec (15.24 km/sec) 

at an altitude of 190 000 ft (57.91 km). 


TABLE I.- NONGRAY STAGNATION-POINT HEAT FLUX 


Altitude 

Sphere radius 

Heat flux 

feet 

meters 

feet 

meters 

Btu/ft^-sec 

kW/m 2 


Velocity, 50 000 ft/sec (15.24 km/sec) 


190 OOO 

57 910 

0.25 

0.0762 

1 806 

20 500 

190 000 

57 910 

1.00 

.3048 

4 840 

54 900 

190 000 

57 910 

5.00 

1.524 

7 530 

85 500 

190 000 

57 910 

15.00 

4.572 

10 630 

120 600 

220 000 

67 060 

5.00 

1.524 

2 114 

24 000 


Velocity, 41 000 ft/sec (12.50 km/sec) 


140 000 

42 670 

5.00 

1.524 

25 740 

292 000 

220 000 

67 060 

5.00 

1.524 

562 

6 380 


Velocity, 34 000 ft/sec (10.36 km/sec) 


140 000 

42 670 

0.25 

0.0762 

410 

4 650 

140 000 

42 670 

5.00 

1.524 

2 760 

31 300 

190 000 

57 910 

5.00 

1.524 

243 

2 760 
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CONCLUSIONS 


On the basis of the results presented, the following conclusions can be stated: 

1. Converged solutions for the radiative heat flux which impinges on the surface of 
an entry vehicle can be obtained with a first-order time -dependent finite-difference 
method if the grid is refined sufficiently. 

2. Only minor differences occur in the results for the nondimen sional surface heat- 
flux distributions for the gray and nongray absorption coefficient models. 

3. The shape of the nondimensional heat -flux distributions is a strong function of 
the flight velocity and a weak function of the size of the body and the altitude. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., February 19, 1969, 

129-01-03-11-23. 
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THE ENERGY EQUATION 

Equations (1) to (4) can be combined so that the energy equation can be written as 


P 


dh 

dt 


dp 

dt 


= Q 


If the enthalpy h is replaced by the right-hand side of equation (15) and e is assumed 
to be a function of p and p, this equation becomes 


9e' 

9P, 


^E + 

p(^\ 

p 

)p dt 

w p 

p 


dt ^ 


From the first and second laws of thermodynamics, 


de = — dp + T dS 
P 2 

where S is the entropy. From equation (A2), 


(Al) 


(A2) 



(A3) 


If e is assumed to be of the form e = e(p,p(p,S)), the derivative (8e/8p)g can be 
written as 


(d e\ _ /8e\ 

Ws ~ W P + 



(A4) 


Equating the right-hand sides of equations (A3) and (A4) and noting that 


yields 



= a 2 



P _ 


P 


pa^ 



P 


(A5) 
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By using equation (A5), equation (Al) can be written as 


dp a 2 »£ Q_ 

dt dt 


/9e\ 
Pr 
\® Ph 


(A6) 


The derivatives (9e/9p )p and (9e/9p)p can be determined in terms of p, e, 
and y(p,e) by replacing e by the expression 


e = 


P 


r|p,e(p,p)] - 1 p 
and differentiating. The derivative (9e/9p)p is written as 


9e\ 

9pjp 


y - 1 + 


9y 


19 log e ei 


and the derivative (9 e/9p)„ is 


9e 

9p 


\ lo g e e /pj 

With the aid of equation (A7) , the energy equation (A6) can be written as 


y - 1 + 


y - 1 + 


9 y 


0 log e 9} 


9 y 


(A7) 


(A8) 


dp 2 dp _ 
dF " a dt ~ 



dy 


v 9 log e e i 


Q 


Equations (A5), (A7), and (A8) can be used to write an expression for the square of the 
speed of sound as 


a 


2 = e 



1) 


y + 


!£ \ 

9 lo ge e Jp 


+ 


dy 


V 9 l°gp PJ 


] 
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CURVE FITS FOR THE THERMODYNAMIC FUNCTIONS 
y(p,e) AND T(p,p) 

The data of Allison (ref. 17) and Browne (ref. 18) have been used to construct curve 
fits for the thermodynamic functions y(p,e) = k and T(p,p). 

The function y = s for equilibrium air for 10 g £- g 10 and - = 1500 is 

e ^o ■ KA o 

given by the following equations: 

y = 1.405 

y = 1.5055 - 0.1255 log 10 ^ 


3.255 - 2.278 log 10 -gA 

^ (B2) 

1 - 0.822 log 1Q J 


0.801 < log 10 = 2.300; 


lQ gio 


(•°S10 HTJ s °' 801 ) < B1 > 


y = ^1.6370 - 0.0404 log 1Q - ^0.2175 - 0.0332 log 1Q £-jlog 10 -jL- 
^0. 1366 - 0.0366 log 10 Aj - ^ 0.0833 - 0.0248 log 10 ^log 10 


1 + exp 


-18.30^1og 10 - 1.943 - 0.032 log 10 jgj 


0.801 < log 10 ^ § 2.300; log 10 §- < 


3.255 - 2.278 log 10 


10 p 


o 1 - 0.822 log 1Q ^r 


(B3) 
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y = (1.5004 - 0.0038 log 10 - (0.1342 - 0.0084 log 10 ^JloglO r|^ 


0.3274 + 0.0091 log 10 - (0.1342 + 0.0016 log 10 ^-jlog 10 


1 + exp 


■15. 


• 85 ( log io ~ 2 * 708 “ °- 0320 log 10 ^ 


2.300 < log 10 = 3.176) (B4) 


The temperature (T in °K) for equilibrium air for 6.00 x 10 -2 = § 2.00 x 10 

and lO"^ < < 10”^ is given by the following equations: 

P 0 


V 0.402 logjP - Ro 

log 10 T = 4.211 + 0.068Ro - 0.392 log 10 ^- - Ro v ° 

' ° 1 - exp 


50.25(log 10 ~~ - R 0 ] 


lo Sl0 7T * R o ( B5 ) 


log 10 T = 4.203 + 0.068 log 10 


{log 1Q j£- = Ro) (B6) 


where 


Ro = 0.941 log 10 ^ 1. 

F o 


970 


(B7) 
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THE PERFECT GAS MODEL 

i_ 

It is assumed that the ratio y = — is constant. By use of this assumption, the 

G 


energy equation 


T dS = dh - ~ dp 
P 


and equations (18) can be combined to obtain 


(y - l) T dS = i dp - y £- dp 


(Cl) 


If the equation for the temperature is of the form 



where the exponents R and S are different from one, equation (Cl) can be written as 


(y - 1)C dS= — ^5. — i- d(p 1- ^) + — ^p 1_ ^d/-J— \ (C2) 

1-Rpl-S 1-S [ p usj 


Let subscripts 1, 2, 3, and 4 denote the properties at the points shown in figure 22. 



Figure 22.- Schematic of dependence of entropy on pressure and density. 
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When equation (C2) is integrated from the states p^, to pg, pg along the 
paths 1,2,3 and 1,4,3 as shown in the figure, the following results are obtained: 


and 


(y - X)C(S 3 - S t ) = - p}-“) 

Pi 



(C3a) 


(r-i)c(Sg 



1 - R 1-S 
P 4 


1-R 1-R 

p 3 p 4 


(C3b) 


Noting that pg = p^, Pg = Pg, P 4 = pg, and p^ = p^, one can equate the right-hand sides 
of equations (C3) and obtain 


1-R 
P _1 

1-S 

^ p l 


p 2~ R , p 3~ R 



(C4) 


Since the first factor in equation (C4) does not vanish in general, the following relationship 

iv r«v 

must exist between y, R, and S: 


1 - S = y(l - r) 

From equation (Cl), the expression for the speed of sound for the perfect gas model 
can be obtained. This equation can be rewritten as 

dp = y jijdp + (y - l) pT dS 
The equation for the square of the speed of sound is 



APPENDIX D 


THE RADIATION INTENSITY, HEAT FLUX, AND NET HEAT INPUT PER UNIT 
VOLUME FOR AN INFINITE, ONE -DIMENSIONAL SLAB OF GAS 

Radiation Intensity 

The differential equation for the radiation intensity l v for a one -dimensional slab 


- £224 - BM (Di) 

where y is the coordinate measured along a line normal to the slab, and cp is the 
angle between the line of propagation of radiation and the y-axis. This coordinate system 
is shown in figure 23. The y coordinate increases from the body to the shock, and the 
y coordinates of the body and the shock are y^ and y s , respectively. 



Figure 23.- Geometry of one-dimensional radiating slab. 

The integration factors expjsec ^Aj^y^yjj and expjsec(i 7 - 0)At-^y,y s ^j are used 

to integrate equation (Dl) for 0 ~ <p = v/2 and v/2 = (p = n, respectively, where the 
integral Aj/(r,s) is given by equation (40). If it is assumed that the body surface does 
not reflect, the solution to equation (Dl) can be written as 
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i£(y,0) = i^(yb) ex P ~ sec <Mi/( y b *y) 


L 

pS=y 

sec 0 \ Q!^(s)Bi/(s)exp[^sec 0Aj,(s,y£]ds 

J s=y b 


(0 i <M 7r/2) (D2a) 


t^(y,0) = I^y s )exp|^sec(7r - 0)A^y,y s jj 

p s =y s 

+ sec(77 - 0) \ o? i; (s)B^(s)exp[^sec(7T - 0) Aj,(y,s)]ds (7 t/2 ^ 0 ^ tt) (D2b) 

dc— 17 


where I^(y^) and I^y s ) are intensities of the radiation entering the slab from the 
body and the gas in the free stream, respectively. 

Heat Flux 

In general, the monochromatic heat flux across an elemental area is given by the 
equation 


iv = -L 


I v cos 


<p do 


(D3) 


If an infinite slab is being considered, the heat flux across an elemental area which is 
parallel to the sides of the slab is given by the equation 


q*»(y) = 




)ify=0 


n(t>=Tl/2 n(p=TT 

\ lr/(y,0)cos cp sin <p dcp + \ I l/ (y,</!))cos <p sin <p dcp 

J 0= 0 J cp=ir/2 


d xp (D4) 


where \p is the azimuthal angle. Since l£ and 1^, are independent of the angle i//, 
the integration with respect to that variable is straightforward. 

It can be seen from equations (D2) that ij and 1^ are functions of sec <p and 
sec( 7 r - cp), respectively. When sec cp and sec(7r - cp) in the first and second integrals 
with respect to cp in equation (D4), respectively, are replaced by w, that equation can 
be written in the form 


Qj/y) = 2 xr 


n CO=oo n (jL)=OC 

f 4(y,^ - e 

^00=1 


ijXy,«) 
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or 


qJy) = 2 tt 


fr b ) r " .-MMs§ . i-(y s )f 7 .-Mw.)S| 

\ ) ^00=1 C O'* ^C0=1 (j)'* 


ns=y ( 

\ a^(s)Bj,(s) l 

J s=y b V J co=l 


a>=°° 


-o)A^(s,y) dcu 


Ids 




r s=y s / pu>=« 

- \ S «„(s)Bj,(s) \ 
J s=y \ J co=l 


w -°° e -o)Aj,(y,s)dco> 
2; 


ds 




(D5) 


By using equation (41) to replace the integrals with respect to on with exponential inte- 
grals, equation (D5) can be written as 


q^y) = 27r^y b )E 3 A„(y b ,y)J - I^y s )E 3 jAj,(y,y s )j + J g _ y B ^( S ) E 2 


A^(s,y) 


“j/(s)ds 


ps=y 

Bj,(s)E 2 

J s=y 


A^y,s) 


o' J/ (s)ds) 


It can be shown from equations (39) and (40) that 


dE n f “1 

^tL a H 

= E^Ja^s^o^s) 

dE n r "I 

d.M 

= “E n _ i ^^(y ,s)J ot^s) 




(D6) 


By using equations (D6), the equation for the monochromatic heat flux is 
q^y) = 2ir/l+(y b ) E 3 Aj,(y b ,y) - C(y s ) E 3 |A^y,y s j 


+ t Y Bj,(s)d /E 3 [A i ,(s,y)] > l + C Ys Bj,(s)d<E 3 

J s=y b l L J J s=y L 


jJa^s) 


(D7) 


The total radiative heat flux can be obtained by integrating the monochromatic heat 
flux over all frequencies: 
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s*V=°° 

q(y) = \ q^y^ 

J v=0 


(D8) 


Net Heat Input to a Volume Element 

The net heat input to a volume element due to monochromatic radiation is equal to 
the negative of the divergence of the monochromatic heat-flux vector -V • 


Qj, = -V • q„ = a v / j I„ dfi - 4 ttB„ 


477 

\ / 

For an infinite slab, Q^(y) is related to q^(y) by the equation 

dq^y) 


(D9) 


Qi/(y) = -• 


dy 


(DIO) 


The expression which is used in this report to determine the net heat input can be 
obtained either from the right-hand side of equation (D9) or from equation (DIO). The 
latter equation is used here. 

When equation (D7) is integrated by parts, the following expression for q^y) is 
obtained: 


q»/(y) = M I t(yb) - B ^(yb) E 3 A ^(yb>y) - ^(y s ) - B ^(ys) E 3 A ^(y>y s ) 


ps=y dB„(s) 
j ds ^3 

-Ms,y) 

p s =y s dB^(s) 
ds " J ds E 3 

A i/y,s) 

J s=y b 

— 


L J 


ds'; 


(Dll) 


From equations (39) and (40), it can be shown that 

dE n 


dy L 
dE n 


A ^(s,y) 


= -E n -l 


dy L 


A ^(y,s) 


= E n-1 


Aj/(s,y) 

A iXy,s) 


My) 

My) 




(D12) 


By using equations (D12) and Leibnitz's rule for differentiating an integral, equation (Dll) 
can be differentiated with respect to y to obtain 
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dq^y) 

dy 


= - 2 Wy)^(yb) - B ^(yb)j E 2^(y b >y)J + ^(y s ) - B ^(y s )j E 2^(y,y s )] 

r s=y dBjy(s) f ^ ps=y g dB„(s) f ^ 

-J s =y h -te- E 2|A^. y )) d = + } -JE- E 2 |A„(y,s)j 


(D13) 


By substituting the right-hand side of equation (D13) into equation (DIO) and using 
equations (D6) on the resulting expression, the following expression for the net monochro- 
matic heat input for a volume element in an infinite slab is obtained: 


Q„(y) = 27rov(y) 


y 

S=Vi_ O' 


^(y b ) - B ^(y b )j E 2^(y b ,y| + j^(y s ) - B E 2 |Aiy(y,y s )j 

EsfA^y^jjj 


s= y i dB j,(s) 


s=y b a v( s ) ds 


d^A.Xs.yj]}-^ 3 - 


s=y s j dB„(s) 


,(s) ds 


(D 14) 


The net heat input due to radiative transfer is obtained by integrating Qj/y) over 
all frequencies: 


Q(y) = C Q„(y)di> 
J v=0 


(D15) 


The monochromatic heat flux can be expressed in a form other than equation (D14) 
by integrating equation (D13) once by parts and substituting the result into equation (DIO). 
The resulting equation is 



(DIO) 
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STAGNATION STREAMLINE SOLUTION FOR GRAY AND 
NONGRAY RADIATION MODELS 

In this analysis, it is assumed that the pressure is constant along the stagnation 
streamline. By using this assumption, the energy equation along the stagnation stream- 
line can be written as 



(El) 


where y is the length coordinate measured along the stagnation streamline. 

Goulard T s stretched coordinate x (ref. 1), which is defined by the differential 
expression 


dx = 



(E2) 


is used in this treatment. As in reference 1, the mass flux normal to the surface is 
approximated by the equation 


P V 

pu = - -—j - — — x (E3) 

5 A 

where and are the density and magnitude of the velocity in the free stream and 

5^ is the shock detachment distance for the adiabatic case. Equations (El), (E2), and 
(E3) can be combined to yield 


-Pqo V qo /P \ 1/2 :: dh _ 

6 A y °s / ^ 



(E4) 


It is assumed that the absorption coefficient can be approximated with a step func- 
tion of the form given in equation (41). Therefore, the divergence of the heat flux is given 
by the expression: 


r 

u — i 


p—oo 

v=0 


M 

Qv dv = ^ Qjq 
N=1 


(E5) 


It is further assumed that the first two steps are optically thin and the remainder are 
optically thick so that 
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«1,«2 < 1 


« N > 1 


(N = 3,4, . . . ,M) 


For the optically thick steps, the expression for the divergence of the heat flux, given by 
equation (D14) or (D16), is simplified with the Rosselind approximation (ref. 11) to yield 


Qn = 


4n 


B 


N,s 0 


36 a (<*n,s 6 a) 9y 


^ a N,s\ 0 ( B N 


°n J 9y \®N,s, 


(N = 3,4, 


,M) (E6) 


where 




(E7) 


and 


B N =r^ 


N+l 


B v dv 


(E8) 


It has been explained previously that the quantities I^y^) and l”^y g j in equations (D14) 
and (D16) are assigned values of zero. In order to approximate equation (D16) for the 
optically thin steps, it is assumed that B^(s) = Bj,(y) in the integral terms of that equa- 
tion. The result is written as 


Qn --^(“N^^BNjS 



(N = 1,2) 


where 


A N [r,s] = A^r,s] = u = "N+l) ( e9 ) 

The exponential integrals in this equation are approximated by the simple exponential 
functions 


-V3z 

~ e 

The exponential functions, in turn, are approximated by the first two terms of their Taylor 
expansions. Therefore, the expressions for the divergence of the heat flux for the opti- 
cally thin steps are written as 
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QN = ^B N>saN , s e A ^L^L 
6 A ’ “N,s b N,s 


- x A N( yb ’ Ys )| 


(N = 1,2) (E10) 


Equation (40) is used to obtain the quantity A N ^y b ,y s ^ in equations (E10). In this treat- 
ment, the contributions to the divergence of the heat flux from the optically thick steps 
Q 3 , Q 4 , . . .,Q^j are neglected. This procedure can be justified if the factors 

B N,s/(°iN s 6 a) and b N,s^“n s 6 a) in equations (E6) and (E10), respectively, satisfy the 
following inequalities: 


B 1,s( Q! 1,s 6 a)> B 2,s(«2,s 6 a) 



(N = 3,4, . . . ,M) 


By using equations (E5), (E10), and the preceding assumption, equation (E4) can be written 
as 


V~h<J— 


1/2 


'oo V OO 


A(A\ = 

dxlhj 


4w 2 b N, S (“n, s «a)^ 
N=1 ' N > £ 


B 


N 


! N,s b N,s 


1 - ^ A N (y b ,y s )J 


(Ell) 


In this treatment, it is assumed that y and p are constant along the stagnation 
streamline. Therefore, it is seen from equations (18) that the density can be expressed in 
terms of the enthalpy as 



(E12) 


An expression for the temperature in terms of the enthalpy can be obtained by combining 
equations (37), (E12), and the assumption of constant pressure to yield 


T_ 

T s 



(E13) 


It is assumed that the absorption coefficients can be approximated by power-law expres- 
sions of the form 


jv = /P_f N m e N 
®N,s \ P sj \ T sJ 


(N = 1,2) (E14) 


Equations (E12), (E13), and (E14) can be combined to yield 
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_"N_ = /h_\ g NS-a N 
“N,s / 


It is shown in figure 24 that the function 
mated by the expression 



given by equation (43) 



for the range 1.15 = log 1r . — = 0.85 where f = 2.42 x 10^ and g = 

1U kT 

tions (E8) and (42) can be used to write and B2 as 




Figure 24.- Approximate solution for integral of Planck function. 


(El 5) 

can be approxi- 

(E16) 

7.45. Equa- 

(E17a) 
(El 7b) 
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where 



From equations (E13) and (E17), it follows that 


and 





b 2,s = B l,s?J 



(E18) 


(El 9) 


(E20) 


By using equations (E12), (E15), (E19), and (E20), equation (Ell) can be written as 


5^) = r l _ 1 ' f A l( Yb ’ y s) 


h \(4+bi)s-a 1+ | 


hS 7 


1 + 


a l,s) 


1 - f A 2 ( y b’ y s) 


1 A i(vy s ) 


' h \(b 2 - b 1 +g)§ - a 2 +a x I 


iP- 


(E21) 


where 


r l = 


_ 4 ( a l,s 6 A) 


oTc 


PocVoohs 


(E22) 


If a new dependent variable £ of the form 

/if i+ H 4+b ~ i) § 

? =■ 


[ 4 + w)s -h-\ 1 - ^ A i(y b ’y s ) 


(E23) 
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is introduced in equation (E21) and the roles of £ and the independent variable x are 
reversed, the energy equation can be written as 


where 


dx 1 

d? + l + 4' 


~ _ (*>2 ~ h + g) g ~ a 2 + ^1 


(4 + bjs - 


" 1-2 


and 

(%s\ 

1 - y- A 2(yb>y s ) 

f* 

\ a l ,s) 

c = 

^1 -^ A l(y b ,y s ) 

1 

£ ( 4 + bj)s - a 



(E24) 


(E25) 


(E26) 


Equation (E24) can be integrated by parts to yield 


x = 6 a exp / 


1 

?s 

\N + 1 

1 + 6 

r N 
_ s s 


F 1,1;2 +i; 

N 




1 + -± 

L 


i + 4 


F 1,1;2 +i; 

N 


-N 


-i-iN 


1 + ■ 


-N 


where the hypergeometric series is written as 



(E27) 


(E28) 


Two of the parameters in equation (E26) for C are not known initially. These param- 
eters are the optical depths A^y^yJ and A 2 (y b ,y s ). Therefore, the final solution 
given by equation (E27) must be determined by iteration. Once a preliminary enthalpy 
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distribution has been obtained, the next approximate values for the optical depths are 
determined from the equations: 



S i x=6 A ~ ~ ~ 1 

AA» N S-a N+5 

x=0 ' S ' 


(N-= 1,2) 


The process is repeated until the successive values for the optical depths converge. The 
final distribution h(y) is determined from the equation 


y = y b 



dx 


(E29) 


It can be seen that equation (E27) yields the solution h = 0 at the stagnation point x = 0. 
A nontrivial value for the stagnation-point enthalpy is obtained in this treatment by extrap- 
olating the values of the enthalpy h at x = 0.055^ and x = 0.16^. 

The results for a gray case and a two-step nongray case are compared with those of 
Olstad (ref. 4) in figure 25. 



Figure 25.- Comparison of results of present initial solution and those of reference 4 for 
stagnation-streamline enthalpy distribution. 
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APPROXIMATE SOLUTION FOR SURFACE HEAT FLUX 

This solution is derived for shock layers for which most of the radiation which 
strikes the body surface comes from the optically thin portion of the spectrum with 
X > 1130A. Equation (45) for the radiative heat flux to the body surface is approximated 
by the expression 


q(y b ,0) = " 2ct f S «i( s >©) T 4 (s,©)ds 

v ' J s=y b 


The dependence of the temperature on the pressure and density and the dependence of the 
absorption coefficient on density and temperature are given by equations (37) and (E14), 
respectively. Therefore, the equation for the heat flux can be written as 


q(©) = 


-20a, 6(0) Tg 5® f S 1 L Ps J - ds 

1,8 6(0) J 3= o r P ( 0 ? s)]( 4+B l) g - a l 

_ Ps _ 


(FI) 


The shock-layer thickness is approximated by the expression 

6(0) = 6(0) (cos 0) -d 


(F2) 


where d is a positive number, and the pressure and density are approximated by expres- 
sions linear in Y of the forms 


= Y cos 2 (0 - e) + (1 - Y)cos20 (F3) 


P(Y,0) _ Y + — (1 - Y)(cos 0) 2 /^ (F4) 

p s p s 

The variation of the density along the body surface is determined by assuming that y is 
constant and that the flow along the surface is isentropic. These conditions are fairly 
good approximations for the cases treated. Since the angle e in equation (F3) is small, 
it is assumed that 
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cos e « 1 

Equations (F2) and (22) are used to obtain an approximate expression for sin e : 

d<5 / d ° „ d6(0)sin 0 


(F5) 


sin e ~ tan e = 


r b + 6 


r b + 6 (°) 

By use of equations (F5) and (F6), equation (F3) can be written as 

Eg£U cosset - 2ae l°) Y Sin2e 
p s r b + 6(0) 

In this treatment, the average value of the pressure across the shock layer is used so 
that 


(F6) 


. cos 2 e + Sa (?h sin ^6 « 1 - B 2 + - 1 6 <°> e 2 « (cos 6 


r b + 6(0) 


d5(0) 


r b +6(0) 


r b + 6(0) 


(F7) 


The integration of the linear expression for the density, which is given by equation (F4), 
can be performed exactly. The result is written as 


r 8 " 1 [ei&si 


ai-(4 + bi)< 


ds = 


^(cos ©) 2 /y 


— * 1+a^ - ^4+bjJs 


i . ^(cos e) 2 /? 


1 + a 


!-(4 + bi)s 


(F8) 


By using equations (FI), (F2), (F7), and (F8), the approximate expression for the surface 
heat flux can be written as 


1+a-i 


< 

q(0) _ 

I 1 ' 

2£i(cose) Vy 

Lp s J 

— \ 


E-( 

'P st \ 1+3r i-( 4+S i) § " 

,<V 

i 

1 - ■ 




(cos 0) 


1 d6(0) 

r b +6(0) 


^4+bjjR-d 


(F9) 


The stagnation -point heat flux is given by the expression: 
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^st 

2 “ 1 , s «(0)ot| 



(F10) 


Equation (F9) contains eight parameters. The exponents R, S, aq, and tq, and the 
ratio y - ^ can be obtained from thermodynamic and absorption coefficient data. The 
present results show that the exponent d which governs the variation of the shock-layer 
thickness with cos 0 can be obtained from adiabatic results. The ratios -=■ 

L r 

and P s t/P s can be determined from a stagnation streamline treatment. For cases for 
which the stagnation streamline density profile varies nonlinearly near the body surface, 
a least-squares procedure should be used to obtain the best linear fit for the density 
distribution. 


6 ( 0 ) 

+ 6 ( 0 )] 


The present treatment has not been extended to include the optically thick radiation 
from the part of the spectrum with X < 1130A because most of the radiation in this spec- 
tral region which strikes the body surface is emitted close to the surface. This radiation 
is affected strongly by the density and temperature gradients near the body surface, which 
are not represented well by this treatment. 
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